#Nature's Kidneys: the Role of Wetland Reserve Easements in Restoring Water Quality
#Nicole Karwowski and Marin Skidmore
#Replication File
#2025

#Figures and Maps

#load libraries
library(ggplot2)
library(dplyr)
library(raster)
library(sf)
library(stringr)
library(haven)
library(tidyverse)
library(RColorBrewer)
library(scales)
library(maps)
library(mapview)
library(leaflet)
library(htmlwidgets)
library(tmap)
library(tigris)
library(cartography)
library(lubridate)
library(fixest)
library(plm)
library(readr)
library(stargazer)
library(ggfixest)
library(cowplot)
library(broom)
library(openxlsx)
library(readxl)

#set wd
setwd("C:/Users/n57m645/OneDrive - Montana State University/UW_Madison/Water quality/Replication Package")

#open huc shapefile
huc8<-st_read("Data/WatershedBoundaries/Raw/HUC8_CONUS/HUC8_US.shp")
huc8$huc2<-substring(huc8$HUC8,1,2)

#keep MRB hucs
huc8<-subset(huc8, huc8$huc2 %in% c("05", "06", "07", "08", "10", "11"))

#Get easements by huc8 totals
load("Data/NRCS easements/Processed/easements_panel_huc_ym.RData")
ease_huc12<-subset(ease_huc12, select=c(huc12, YearMonth, RestoredEasementAcres))
ease_huc12$huc8<-substring(ease_huc12$huc12, 1,8)
ease_huc12$year<-substring(ease_huc12$YearMonth, 1, 4)
ease_huc12<-subset(ease_huc12, as.numeric(ease_huc12$year)<2019)
sum(ease_huc12$RestoredEasementAcres, na.rm=TRUE)

ease_huc8 <- ease_huc12 %>%
  group_by(huc8) %>%
  summarise(EasementAcres = sum(RestoredEasementAcres, na.rm = TRUE))
sum(ease_huc8$EasementAcres)

ease_huc8<-distinct(ease_huc8)
colnames(ease_huc8)<-c("huc8", "Easements")
summary(ease_huc8$Easements)

#merge to shapefile
huc8<-merge(huc8, ease_huc8, by.x=c("HUC8"), by.y=c("huc8"), all.x=TRUE, all.y=FALSE)

#if easements are missing, make it 0. 
huc8$Easements<-ifelse(is.na(huc8$Easements), 0, huc8$Easements)
sum(huc8$Easements)

rm(ease_huc12)
rm(ease_huc8)
summary(huc8)

#Get water quality measures
#open panel
load("Data/Main/main_huc12_panel.RData")

#keep only relevant years 
merge<-subset(merge, as.numeric(Year)>=1985)

#Get water quality readings
#keep only observations that have an ammonia or phosphorous measurement
merge<-subset(merge, !is.na(ammonia_filtered) | !is.na(TP_filtered) | !is.na(TKN_filtered))


#pull stats that i need
stats<-merge %>% group_by(huc8) %>% mutate(Ammonia=mean(ammonia_filtered, na.rm=TRUE),
                                           TKN=mean(TKN_filtered, na.rm=TRUE),
                                           Phosphorus=mean(TP_filtered, na.rm=TRUE),
                                           Treated=max(EverTreated, na.rm=TRUE))

stats<-subset(stats, select=c(huc8, Ammonia, TKN, Phosphorus, Treated))
stats<-distinct(stats)
summary(stats)

stats <- stats %>% mutate(across(everything(), ~ ifelse(is.nan(.), NA, .)))

#merge marb hucs and ammonia and phosphorous levels
huc8<-merge(huc8, stats, by.x="HUC8", by.y="huc8", all.x=TRUE, all.y=FALSE)
colnames(huc8)
summary(huc8)

#convert easements to percentage share of huc8. 
huc8$EasementsP=100*huc8$Easements/huc8$AREAACRES
summary(huc8$EasementsP)

sum(huc8$Easements, na.rm=TRUE)
#there are 1,758,137


# Create categorical variable for EasementsP
huc8$EasementsP_cat <- cut(huc8$EasementsP, 
                           breaks = c(-Inf, 0, 0.001, 0.01, 0.1, 1, 8),  # Defined bins
                           labels = c("0", "0 - 0.001", "0.001 - 0.01", "0.01 - 0.10", "0.10 - 1.00", "1.00 - 8.00"),
                           include.lowest = FALSE)
table(huc8$EasementsP_cat)

# Convert to factor for ordering
huc8$EasementsP_cat <- factor(huc8$EasementsP_cat,
                              levels = c("0", "0 - 0.001", "0.001 - 0.01", "0.01 - 0.10", "0.10 - 1.00", "1.00 - 8.00"))
table(huc8$EasementsP_cat)

custom_palette <- c("0" = "#D9C2A3",      
                    "0 - 0.001" = "#E5F5E0", 
                    "0.001 - 0.01" = "#C2E699", 
                    "0.01 - 0.10" = "#78C679", 
                    "0.10 - 1.00" = "#31A354", 
                    "1.00 - 8.00" = "#006400") 

#nutrient color palette
nutrient_palette <- c("#E6F2FF", "#99CCFF", "#3399FF", "#0066CC", "#003366")

#create treated sample subset
#this is where we have wq data and where we have easements. 
treated<-subset(huc8, huc8$Treated==1)
treated1<-st_union(treated)

treated_bold<-tm_shape(treated1) + tm_borders(col = "maroon", lwd = 2, lty="dashed")


# Load the USA boundaries
usa<-states(cb=TRUE, resolution="500k", class="sf")
usa<-st_transform(usa, crs=st_crs(huc8))
usa<-st_crop(usa, st_bbox(huc8))

#USA map
base<-tm_shape(usa)+tm_polygons(col="lightgray")+tm_borders(col = "white", lwd = 1)

#Nutrient levels maps
ammonia_map<-tm_shape(huc8) + tm_polygons(fill="Ammonia", fill.scale = tm_scale_intervals(style="quantile", values = nutrient_palette), fill.legend = tm_legend(title = "Ammonia", position = tm_pos_in("right", "top")))+tm_scalebar()
tkn_map<-tm_shape(huc8) + tm_polygons(fill="TKN", fill.scale = tm_scale_intervals(style="quantile", values = nutrient_palette), fill.legend = tm_legend(title = "TKN", position = tm_pos_in("right", "top")))+tm_scalebar()
phosphorus_map<-tm_shape(huc8) + tm_polygons(fill = "Phosphorus", fill.scale = tm_scale_intervals(style = "quantile", values = nutrient_palette), fill.legend = tm_legend(title = "Phosphorus",  position = tm_pos_in("right", "top")))+tm_scalebar()

#Easement maps
ease_map<-tm_shape(huc8) + tm_polygons(fill="EasementsP_cat", fill.scale = tm_scale_categorical(values=custom_palette), fill.legend = tm_legend(title = "Wetlands (%)", position = tm_pos_in("right", "top")))+tm_scalebar()

#put it all together
m1<-base+ammonia_map+treated_bold
m2<-base+tkn_map+treated_bold
m3<-base+phosphorus_map+treated_bold
m4<-base+ease_map+treated_bold

#create panel at the end
all<-tmap_arrange(m1, m2, m3, m4, ncol=2, nrow=2)

#save it
tmap_save(all, filename="Figures/jaere_maps_nutrients_wrp_revised.eps", width=16, height=12)
tmap_save(all, filename="Figures/jaere_maps_nutrients_wrp_revised.jpg", width=16, height=12)

################################################################################
#Benefits by PWS. 

#Explore the PWS. 
load("Data/Water System EPA/Processed/count_pws_type.RData")

#create count of benefits
fc_huc8$benefits_size<-1289*fc_huc8$n_pws_active_sw_small+2072*fc_huc8$n_pws_active_sw_medium+17265*fc_huc8$n_pws_active_sw_large
benefits_size<-subset(fc_huc8, select=c(huc8, n_pws, benefits_size))

huc8<-merge(huc8, benefits_size, by.x="HUC8", by.y="huc8", all.x=TRUE, all.y=FALSE)
huc8$benefits_size<-ifelse(is.na(huc8$benefits_size), 0, huc8$benefits_size)


huc8$ihs_Easements<-asinh(huc8$Easements)


huc8$benefits_size_all<-huc8$benefits_size*huc8$ihs_Easements

huc8$n_pws<-ifelse(is.na(huc8$n_pws), 0, huc8$n_pws)

summary(huc8$benefits_size)
summary(huc8$benefits_size_all)


#add one wetland
benefits_bypws_one<-base+tm_shape(huc8) + tm_polygons(fill="benefits_size", 
                                          fill.scale = tm_scale_intervals(style="fixed", 
                                          breaks = c(0, 1, 1000, 10000, 100000, Inf),
                                          labels=c("0", "0-1,000", "1,000-10,000", "10,000-100,000", "100,000-2,087,000"),
                                          values = c("#FFFFCC", "#FFDD66", "#FFAA00", "#3366CC", "#003366")),
                                          fill.legend = tm_legend(title = "Potential Benefits {$}", na.label="0", position = tm_pos_in("right", "top")))+tm_scalebar()


#save figures
tmap_save(benefits_bypws_one, filename="Figures/jaere_maps_benefits_bypws_one_low.eps", width=8, height=6)
tmap_save(benefits_bypws_one, filename="Figures/jaere_maps_benefits_bypws_one_low.jpg", width=8, height=6)


#benefits from all wetlands. 

benefits_bypws_all<-base+tm_shape(huc8) + tm_polygons(fill="benefits_size_all", 
                                                      fill.scale = tm_scale_intervals(style="fixed", 
                                                  breaks = c(0, 1, 1000, 10000, 100000, 1000000, Inf),
                                                  labels=c("0", "0-1,000", "1,000-10,000", "10,000-100,000", "100,000-1,000,000", "1,000,000-18,564,000"),
                                                  values = c("#FFFFCC", "#FFDD66", "#FFAA00", "#3366CC", "#003366", "#2E1A47")),
                                                  fill.legend = tm_legend(title = "Realized Benefits {$}", na.label="0", position = tm_pos_in("right", "top")))+tm_scalebar()

tmap_save(benefits_bypws_all, filename="Figures/jaere_maps_benefits_bypws_all.eps", width=8, height=6)
tmap_save(benefits_bypws_all, filename="Figures/jaere_maps_benefits_bypws_all.jpg", width=8, height=6)


################################################################################
#Create a map of the benefits gained
#by nutrient elasticity

#for each huc8, create a measure of whether the huc8 is considered low, medium, or high ammonia. 
#use terciles
# Compute tercile breakpoints
tercile_breaks <- quantile(huc8$Ammonia, probs = seq(0, 1, length.out = 4), na.rm = TRUE)

# Assign tercile groups
huc8$Ammonia_Tercile <- cut(huc8$Ammonia, 
                            breaks = tercile_breaks, 
                            labels = c("Low", "Medium", "High"), 
                            include.lowest = TRUE)


#multiply the ihs(Esements) by the WRP VC $ from our cost benefit tables
# for low nitrogen areas, $242.
#for medium nitrogen areas, $1,381
#for high nitrogen areas, $17,265

#benefits by nutrient if we add one wetland easement. 
huc8$benefits_nutrient<-ifelse(huc8$Ammonia_Tercile=="Low", 242, NA)
huc8$benefits_nutrient<-ifelse(huc8$Ammonia_Tercile=="Medium", 1381, huc8$benefits_nutrient)
huc8$benefits_nutrient<-ifelse(huc8$Ammonia_Tercile=="High", 17265, huc8$benefits_nutrient)

#make it 0 for places with 0 PWS
summary(huc8$n_pws)
summary(huc8$benefits_nutrient)
huc8$benefits_nutrient<-ifelse(huc8$n_pws==0, 0, huc8$benefits_nutrient )

#benefits by nutrient times all easements in 2018. 
huc8$benefits_nutrient_all<-huc8$benefits_nutrient*huc8$ihs_Easements

#graph the benefits
summary(huc8$benefits_nutrient)
table(huc8$benefits_nutrient)
summary(huc8$benefits_nutrient_all)

#should only have three categories. since its low, medium and high. 
benefits_byn_one<-base+tm_shape(huc8) + tm_polygons(fill="benefits_nutrient", 
                                                    fill.scale = tm_scale_intervals(style="fixed", 
                                                                                    breaks = c(0, 242, 1381, 17265, Inf),
                                                                                    labels=c("0", "242", "1,381", "17,265"),
                                                                                    values = c("#FFFFCC","#FFDD66", "#FFAA00", "#3366CC")),
                                                    fill.legend = tm_legend(title = "Potential Benefits {$}", position = tm_pos_in("right", "top")))+tm_scalebar()

tmap_save(benefits_byn_one, filename="Figures/jaere_maps_benefits_byn_one.eps", width=8, height=6)
tmap_save(benefits_byn_one, filename="Figures/jaere_maps_benefits_byn_one.jpg", width=8, height=6)


benefits_byn_all<-base+tm_shape(huc8) + tm_polygons(fill="benefits_nutrient_all", 
                                                    fill.scale = tm_scale_intervals(style="fixed", 
                                                                                    breaks = c(0, 1, 1000, 10000, 100000, Inf),
                                                                                    labels=c("0", "0-1,000", "1,000-10,000", "10,000-100,000", "100,000-204,100"),
                                                                                    values = c("#FFFFCC", "#FFDD66", "#FFAA00", "#3366CC", "#003366")),
                                                    fill.legend = tm_legend(title = "Realized Benefits {$}", position = tm_pos_in("right", "top")))+tm_scalebar()

tmap_save(benefits_byn_all, filename="Figures/jaere_maps_benefits_byn_all.eps", width=8, height=6)
tmap_save(benefits_byn_all, filename="Figures/jaere_maps_benefits_byn_all.jpg", width=8, height=6)

################################################################################

#what is the total sum of benefits? using the PWS and all wetlands
sum(huc8$benefits_size_all)
#194,261,825

################################################################################
#Plot payback period

# Functions
small <- function(x){335000-1289*x}
medium <- function(x){335000-2072*x}
large <- function(x){335000-17258*x}
highN <- function(x){335000-17265*x}
medN <- function(x){335000-1381*x}
lowN <- function(x){335000-242*x}

# Theme
My_Theme <- theme(
  title = element_text(size = 12),
  text = element_text(color = "black"),
  strip.text = element_text(color = "black", size = 14),
  axis.title.x = element_text(size = 14),
  axis.text.x = element_text(size = 14),
  axis.text.y = element_text(size = 14),
  axis.title.y = element_text(size = 14)
)

# Create data frame
x_vals <- seq(0, 100, length.out = 200)
plot_df <- data.frame(
  x = rep(x_vals, 6),
  y = c(lowN(x_vals), small(x_vals), medN(x_vals),
        medium(x_vals), large(x_vals), highN(x_vals)),
  group = factor(rep(c("Low nitrogen", "Small facility", "Medium nitrogen",
                       "Medium facility", "Large facility", "High nitrogen"), each = length(x_vals)),
                 levels = c("Low nitrogen", "Small facility", "Medium nitrogen",
                            "Medium facility", "Large facility", "High nitrogen"))
)

# Add linetype and shape
plot_df <- plot_df %>%
  mutate(
    linetype = case_when(
      group == "Medium facility" ~ "dotted",
      group == "High nitrogen" ~ "dashed",
      TRUE ~ "solid"
    ),
    shape = case_when(
      group == "Low nitrogen" ~ 16,
      group == "Small facility" ~ 15,
      group == "Medium nitrogen" ~ 17,
      TRUE ~ NA_real_
    )
  )

# Define colors
colors <- c(
  "Low nitrogen"="#E9C46A", 
  "Small facility"="#F4A261", 
  "Medium nitrogen"="#d7642c", 
  "Medium facility"="#2A9D8F", 
  "Large facility" ="#00a0e1", 
  "High nitrogen"="#466eb4"
)

# Subset points for markers to space them out
marker_df <- plot_df %>% filter(row_number() %% 10 == 0)

# Plot
plot <- ggplot(plot_df, aes(x = x, y = y, group = group, color = group)) +
  geom_line(aes(linetype = linetype), size = 1.5) +
  geom_point(data = marker_df, aes(shape = shape), size = 3, na.rm = TRUE) +
  scale_color_manual(name = "Legend", values = colors) +
  scale_linetype_identity() +
  scale_shape_identity() +
  guides(color = guide_legend(
    override.aes = list(
      linetype = c("solid", "solid", "solid", "dotted", "solid", "dashed"),
      shape = c(16, 15, 17, NA, NA, NA)
    )
  )) +
  xlab("Payback period (years)") +
  ylab("Wetland restoration investment ($)") +
  scale_x_continuous(limits = c(0, 100), breaks = seq(0, 100, 10)) +
  scale_y_continuous(labels = comma, limits = c(-1000, 350000)) +
  theme_bw() + My_Theme

plot

ggsave(plot, filename="Figures/payback_graph.eps", width=8, height=5, units="in")
ggsave(plot, filename="Figures/payback_graph.jpg", width=8, height=5, units="in")

#################################################################################
#Easement acreage over time
load("Data/NRCS easements/Processed/easements_panel_huc_ym.RData")

#keep only MRB hucs
ease_huc12$huc2<-substring(ease_huc12$huc12, 1, 2)
ease_mrb<-subset(ease_huc12, ease_huc12$huc2 %in% c("05", "06", "07", "08", "10", "11"))
ease_mrb$Year<-substring(ease_mrb$YearMonth, 1, 4)
ease_mrb<-subset(ease_mrb, select=c(huc12, Year, YearMonth, RestoredEasementAcres))

#create MRB by year panel
ease_mrb_year <- ease_mrb %>%
  group_by(Year) %>%
  summarise(
    total_restored = sum(RestoredEasementAcres, na.rm = TRUE)
  ) %>%
  arrange(Year) %>%
  mutate(cum_restored = cumsum(total_restored))

#plot
library(scales)

easement <- ggplot(ease_mrb_year, aes(x = as.numeric(Year), y = cum_restored)) +
  geom_area(fill = "darkgreen") +
  geom_line(color = "black", size = 1.2) +
  geom_point(color = "black", size = 2) +
  theme_bw(base_size = 16) +
  theme(
    panel.grid.major = element_line(color = "grey85"),
    panel.grid.minor = element_blank(),
    axis.text = element_text(color = "black"),
    axis.title = element_text(face = "bold"),
    plot.title = element_text(face = "bold", hjust = 0.5)
  ) +
  labs(
    x = "Year",
    y = "Cumulative Restored Wetlands Acres"
  ) +
  scale_y_continuous(labels = comma)

easement
ggsave(easement, file="Figures/cumulative_wetland_acres.jpg", width=8, height=5, units="in")
ggsave(easement, file="Figures/cumulative_wetland_acres.eps", width=8, height=5, units="in")

